static void InitSRBAngCoords ()
{
  VecR e;
  real eAng[ 3 ], phi, theta, psi, cphi, ctheta, cpsi, sphi, stheta, spsi;
  int n;

  DO_SITE {
    if ( site[ n ].mechType == SRB ) {
      VRand( &e );
      eAng[ 0 ] = atan2( e.x, e.y );
      eAng[ 1 ] = acos( e.z );
      eAng[ 2 ] = 2. * PI * RandR();
      /* 3/9/05 - assign rotation matrix for symplectic-splitting DLM integration - from allen-tildesley p 86 */
      phi   = eAng[ 0 ];
      theta = eAng[ 1 ];
      psi   = eAng[ 2 ]; /* THIS IS WORKING, ROTATION MATRIX IS (ORTHONORMAL AND) EQUAL TO THE ONE OBTAINABLE BY THE QUATERNIONS */
      cphi = cos( phi );
      sphi = sin( phi );
      ctheta = cos( theta );
      stheta = sin( theta );
      cpsi = cos( psi );
      spsi = sin( psi );
      siteRotationMatrix[ n ].u[ 0 ] =  cphi * cpsi - sphi * ctheta * spsi;
      siteRotationMatrix[ n ].u[ 1 ] = -cphi * spsi - sphi * ctheta * cpsi;
      siteRotationMatrix[ n ].u[ 2 ] =  sphi * stheta;
      siteRotationMatrix[ n ].u[ 3 ] =  sphi * cpsi + cphi * ctheta * spsi;
      siteRotationMatrix[ n ].u[ 4 ] = -sphi * spsi + cphi * ctheta * cpsi;
      siteRotationMatrix[ n ].u[ 5 ] = -cphi * stheta;
      siteRotationMatrix[ n ].u[ 6 ] =  stheta * spsi;
      siteRotationMatrix[ n ].u[ 7 ] =  stheta * cpsi;
      siteRotationMatrix[ n ].u[ 8 ] =  ctheta;

      GetOrientPrincAxisZ_FromRotMat( &siteOrientationVector[ n ], site[ n ].qMat ); /* convert rotation matrix to orientation vector */

    }
  }
}

